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We discuss the stability of a homogeneous two-dimensional Bose gas at finite temperature against 
formation of isolated vortices. We consider a patch of several healing lengths in size and compute 
its free energy using the Euclidean formalism. Since we deal with an open system, which is able to 
exchange particles and angular momentum with the rest of the condensate, we use the symmetry- 
breaking (as opposed to the particle number conserving) formalism, and include configurations 
with all values of angular momenta in the partition function. At finite temperature, there appear 
sphaleron configurations associated to isolated vortices. The contribution from these configurations 
to the free energy is computed in the dilute gas approximation. We show that the Euclidean action 
of linearized perturbations of a vortex is not positive definite. As a consequence the free energy of 
the 2D Bose gas acquires an imaginary part. This signals the instability of the gas. This instability 
may be identified with the Berezinskii, Kosterlitz and Thouless (BKT) transition. 
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I. INTRODUCTION 



Below a certain temperature in a three-dimensional Bosonic system, long range order is established in 
an ordered phase called the Bose-Einstein condensate (BEC). In two-dimensional (2D) Bose gas, there does 
not exist an ordered state since its existence means that the correlation of its fluctuations is logarithmically 
divergent, as proven by Mermin, Wagner ilj, Hohenberg [2] and Coleman [3]. However, it has been shown 
by Berezinskii [J] , Kosterlitz and Thouless [5] (BKT) that there exists a superfluid phase with a quasi-long- 
range order below a certain temperature. The superfluid phase has only the bounded vortex pairs but above 
the BKT temperature single vortices proliferate as this is the more stable configuration [5j. 

Many theoretical studies on BKT transition in 2D Bose gas, including the original theory of BKT, are 
based on equilibrium thermodynamics of many-body systems [5HlO| . There have been studies of 2D vortex 
dynamics, viewed as massive charged particles in relativistic two-dimensional electrodynamical systems [TTJ 
[12] . This analogy has been applied to the study of vortex dynamics of 2D superfluid via a Fokker-Planck 
equation [T3] and using field theoretical approaches [T^HTB] . There were also numerical studies based on the 
Gross-Pitaevskii equation [T71[TH] on the lifetime of spontaneous decay of a pancake-shaped condensate with a 
vortex [IH] . There were also studies around the critical region of BKT transition with Monte Carlo simulations 
pOl [21] with the local density approximation. With the use of projected Gross-Pitaevskii equation (PGPE) 
|22j . the thermal activation of vortex pairs in the presence of a harmonic trap p5H^ and its emergence of 
superfluidity f^Hj were studied, and the various consequences of the improved mean- field Holzmann-Chevalier- 
Krauth (HCK) theory [57] of the 2D Bose gas ^28.-,31] were presented. The non-equilibrium response of a 
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2D Bose gas is less understood. One encounters such a situation when the trap is suddenly turned off, as is 
done in recent experiments [521 - f5i] described below. 

Recently 2D quantum Bose gas has been experimentally realized by the Dalibard group [32l [33] by slicing 
a 3D BEG into pieces of "pancakes" with ID optical lattices, and by the Phillips group [53] through trapping 
the atoms in a 3D harmonic potential with a very large frequencies in one of the directions. In both 
experiments, measurements on the gas are performed some time after the confining potential is abruptly 
turned off. Dalibard group's experiment showed that there are more isolated vortices formed at higher 
temperatures. The Phillips group measured the density profile after 10 ms time of flight , and identified 
different states of the gas. In one regime the gas develops a bimodal distribution with only thermal and quasi- 
condensate components without long range order, as different from a superfluid. For a sufficiently long time 
of flight, they observe a trimodal distribution with thermal, quasi-condcnsate and superfluid components 
indicative of a BKT transition 

In this paper, we compute the free energy of a 2D Bose gas by means of thermal field theory. We consider 
the action in the Madelung representation (in terms of density and phase), and convert it to a Euclidean 
action by a Wick rotation in time and in phase. The system that we study is a patch of a size of several 
healing lengths within the larger 2D gas. Because we are dealing with the homogeneous configuration, we 
put no confining potential, i.e., ^(x) = 0. Since the vortices form at the center of the gas patch [2S] at 
the beginning where the density of the gas is effectively homogeneous, and the vortex core structure is very 
small compared to the size of the gas patch, we expect to reduce the physically more relevant inhomogeneous 
situation to the homogeneous situation discussed here through a local density approximation. Since particle 
number and angular momentum are not conserved for this system, we do not constraint the former (unlike 
in the particle number conserving formalism, see e.g., [35]) and consider configurations with all values of 
angular momentum. In particular, we consider configurations with different numbers of vortices and the 
fluctuations around them. Although these configurations are time-independent, they have finite euclidean 
action as a consequence of the compactification of the euclidean time axis, namely euclidean time is periodic 
with periodicity These time independent configurations with nonzero angular momentum play in our 
problem the same role as the usual sphaleron configurations in electroweak symmetry breaking [36'. The 
contribution from these configurations to the free energy is computed within the dilute gas approximation. 

We find that the Euclidean action for fluctuations around an isolated vortex is not positive definite. In real 
time, this means an instability of the isolated vortex, and we characterize the direction of greatest instability 
in configuration space. In imaginary time the fact that the Euclidean action is not positive definite means 
that the partition function must be defined by an analytic continuation, whereby the free energy becomes 
complex. We calculate the imaginary part of the free energy due to this instability. This is similar to the 
argument of Langer who considered the decay of a metastable state due to classical fiuctuations [ST] [38] , that 
of Coleman who considered the quantum fiuctuations around the spatially-separated instantons |39H41| , and 
that of Affleck who considered the decay of a quantum-statistical metastable state using instantons [1^ 143] . 
We find that the canonical 2D Bose gas is indeed unstable at finite temperature, and the decay rate, which 
is also the rate of vortex nucleation, increases with temperature. For T > TbkTi the gas evolves to a state 
of isolated vortices. 

The paper is organized as follows. In Sec. [lljwe introduce the Gross-Pitaevskii treatment and write it in 
the Madelung representation. In Sec. Ill we obtain the Euclidean action by a Wick rotation and model the 
density profile of the gas with a vortex at the origin. In Sec. |IV[ we introduce the linear perturbation about 
the configuration for each q. In Sec. |V| we outline the formalism of computing the lifetime of the gas and 
obtain the BKT transition temperature. In Sec. | VI[ we use Bohr-Sommerfeld quantization to show that the 
effective energy is complex, indicating the instability of the 2D Bose gas. We end with conclusions in Sec. 

m 



II. MODEL 



The dynamics of a two-dimensional (2D) Bosonic atomic system with a (5-potential inter-atomic interaction 
is described by the action [TT] [TH] 

S= f dtd'^x lih'^''^ - h\ , (1) 
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where ^'(x) and 5'^^(x) are respectively the annihilation and creation operators of an atom at point x. The 
Hamiltonian is 

H ^ —V^^V^ + F[^'^-^], (2) 

and 

i^[p] = (t/(x)-M)p+^5p'- (3) 
where g is the coupling constant due to the (5-potential between the atoms. In the Madelung representation 

* = ^pe'\ (4) 

the density of atoms in the lowest macroscopically occupied state p and the phase 9 are canonical to each 
other, obeying the commutation relation [HI fi5] 



With Q, the action ([T]) is written as 



where 



and 



[p(x),^(x')] = -z<5(x-x'). (5) 



5= / dtcP^ {ne^-H\, (6) 



H=^{yh6f+F,[p], (7) 



F,[p\=F[p\ + ^p{'^pf- (8) 

The length scale that characterizes the local alteration of the gas density healing back to the mean-field 
density is given by the healing length, which is 

2 

Experimentally there is a harmonic trap to prepare the initial patch of Bosc gas in two-dimensions. At the 
time when the trap is turned off, the Bose gas is still highly inhomogeneous. However, in recent experiments, 
the vortex core structure is very small compared to the patch of quasi-two-dimensional Bose gas. Take Phillips 
group's experiment for example. Sodium atom is used and therefore m ^ 3.8 x 10~^^ kg. And cjj^ = 20 Hz 

and (jjz ^ 1 kHz [34]. It is known that p — gpo, and po = ^ [20 and A = \J~^^x being the thermal de 



Broglie wavelength. Near the transition point, T ~ 100 nK [33l [34] , Po ^ 3 x 10^^ mr'^ . By 5 = -^fj (where 
a is the scattering length and Lz is the thickness of the gas) and the fact that for most current experiments 
f; ~ 5^ [IS], 5 9.7 X 10-"*^ J m2. Then - 2.9 x lO'^^ J. Then from jg]), - 2.5 x IQ-i^ m^. The area 

of the gas is given by the circle of the TF radius, A - 7ri?^ - - 1.20 x 10"^ m^. Hence A >> which 

means the vortex structure is very small compared to the size of the gas. Hence, the experimental situation 
can be recovered from our subsequent analysis through local density approximation j21| that locally the gas 
is effectively homogeneous at the center of the trap [3H], which is best described by ^(x) = 0. 

III. EUCLIDEAN ACTION 

To compute the partition function of such a system, we perform a Wick rotation by writing t — —it. To 
preserve the same canonical relation between the density and the phase ([s]) and to keep the density real, the 
phase has to be rotated accordingly by 

X = -ihe. (10) 
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whence exp (i^) becomes exp (— f )• The action in this EucUdean space is given by 



where the Hamiltonian density is 



5= / drd^x \x^+n}, (11) 



n = -^{Vxf + F,[p]. (12) 
Let us introduce the foUowing dimensionless variables 

T ^ -s,r ^ ^y,p^ -n, (13) 
9 

and the Euchdean phase 

X = (14) 



With theses new variables the EucUdean action (11) becomes 



<|-2„(V.,f-„4 + i^l ,15, 



Because '3/ is a single-valued function its value is unchanged upon having the phase added by 27rg, for 
any integer q, to it does not change the value of the field. As a result, for any integer q, 

dl ■ VyC = -2TTiq, (16) 

where the line integral goes around a loop about a point. If the vorticity q is positive (negative) while the 
loop is small enough, there is a vortex (an antivortex) at that point whereas q — indicates there is no 
vortex at that point. But if the loop of the line integral is larger, q is the sum of the vorticities of all vortices 
inside the loop, while vortex and antivortex cancel each other in the integration. The phase may have a 
curl-free part even if there is a vortex. The simplest configuration representing a single vortex at the origin 
has C = —iqf- The Euclidean angular momentum density of the system is given by 

dx 

I = p— = -ihqp, (17) 

which is proportional to q. The fact that the angular momentum commutes with the Hamiltonian and is 
conserved implies the conservation of vorticity in the whole system. 

Assuming there is a vortex at the origin with density profile nq{y), presumed to be rotationally invariant, the 



equation of motion is obtained by putting ^ — —iqf into ( 15 ) 



For q = 0, Uq — 1 exactly. In the general case, it is convenient to introduce an "Euclidean wave function of 



the condensate" by writing = ip^ [HI |47] . It then becomes 



l_d 
ydy 



dy 



y- 



y2 



(19) 



The vortex solution interpolates between the no-vortex profile ipq = 1 for y oo and the trivial solution 
ipq = ioi y 1-^ 0. Eq. (19) may be solved numerically (see [T51I37]). For large y, we may expand ipq in 
inverse powers of y^: 
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V', = 1 



2y4 



O 



(20) 



Likewise for the density: 



= 1 - 



2(7^ Sq 



(21) 



For y << 1, the cubic term in (19 1 can be neglected, and ipq becomes a Bessel function [4F. For our 
purposes, it is enough to keep only the first (linear) term in the Taylor expansion of Tpg. The density profile 
is then quadratic 



ni = 0.08?/^ 



(22) 



We shall adopt the approximation (21 1 for i/ > 2.7 and (22 1 otherwise. The matching point and the constant 
in ( 22 ) are chosen so the approximated density profile is smooth (see Fig. [T]) 




FIG. 1: The density profile for an isolated vortex at the origin, as given by (22 1 for y < 2.7 and (20 1 for y > 2.7 



IV. LINEAR PERTURBATION 



Consider linear perturbations around a configuration of the 2D Bose gas with a vortex at the origin: 

n^ng{l + S),(:^-iq^ + Ci, (23) 

where 6 — 5{y,Lp,s) and Ci — C,i{y,Lp,s) are functions of the radial and azimuthal coordinates y and (p. 
Define the operator 



1 d 



d 



yuq dy V^"'' dy J ' g^p 



1 



(24) 



Note that for q = (i.e., Uq — 1), — V^. Putting the perturbation (23) in the action (15), it becomes 



siq)^Mq) 



2 /-M/S 



n„ 



1 ~. 



ds I d^ynq5[ f--Vl]5 



lo 

+ ^£'d.sJd^yCi{2nqVl)Ci 
hi,e r\ f,2 , f d ^^q d 



(25) 
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where 



Mq) = ^ r ds fd 

9 Jo J 



2nq5^ 1 f dn,^^ 



2 2n„ \ 



dy ■ yug, 



(26) 



which is the equihbrium free energy, with the second equahty owing to (18). Then the equations of motion 
are given by 



— = 71,(5 + — V (5, 

OS y' dip " 

85 Aiq dS ^y2^ 

ds y"^ dip ^ 

The Fourier transform of the fluctuations can be defined as 

oo 

Siy,<f,s) = Sj{y,s 
i=-oo 

oo 

Ci(y,<yf,s) = Y (ij{y,s)e 

If S and are real, then 



j = -oo 



tjlp 



The representation ( 30 1 assumes that 



5; = <5_„cr, = ci-,. 



dl ■ VCi = 0, 



(27) 
(28) 

(29) 
(30) 

(31) 
(32) 



which means the fluctuation does not change the total vorticity. The counterpart of — in the Fourier 
representation is 



— — (yn — V- 
yuq dy y''dyj y 



(33) 



The action ( 25 1 becomes 



ds J dy yriq ^-j y + l^L^qj 
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7 + ) 



2 /-M/S 



j = -oo 
oo 



/•MP f ^ 



(34) 



2 /-P/S 



j = -oo 
oo 



j ds j dy- yuq ^ Ci- 



3 = -oo 



d_ ^ Aqj 
ds y"^ 



2 / "J- 



This action can be further simplified. Define 



Suppose F and / are related by F = then Ljq and are related by 



(35) 



(36) 
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Define the covariant differential operator, 



as 



(37) 



which can be seen as the time-derivative in a frame corotating with the vortex. With the transformation of 
the fluctuations, 



the action ( 34 ) is then rewritten as 



S{q) sa J^f^(q) + ^^f^ ds J dy ■ y ^ 



(38) 



(39) 



j=-oo 

^-j 2 '"^^i ~ Cl-j(2ijg)Clj + Cl-jDs.JqSj 



From ( 27 1 and ( 28 ) , or from the action ( 39 ) , the equations of motion in terms of the new operators are 



(riq + Ljq)5j — Dg^^jqC,ij, 



(40) 
(41) 



LIFETIME OF THE CONDENSATE 



Consider a Bose gas confined to a region of size L. We define as our system a part of the Bose gas with 
hnear size / smaUer than L but greater than the heahng length ^, i.e., L ^ I ^ The Bose gas within 
this system is interacting with other atoms outside, which act as a reservoir of energy, particle number and 
angular momentum. Therefore the total vorticity q of our system is not conserved. The equilibrium state is 
described by the partition function |49j 



Z= f2 J dSdCi J D[S]D[C,]exp , 



(42) 



g— — oo 



where the periodic boundary conditions [50 S — (5(0) ~ S{fi(3) and Ci = Ci(0) = Ci(/^/5) have been incorpo- 
rated in the evaluation of the path integral. 

Setting an upper cutoff at y = A, the Euclidean action (essentially the free energy divided by ksT) of the 
system with no vortex is obtained by putting no = 1 in ([26]) 



(43) 
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and that with one vortex of vorticity q is obtained after putting the asympotic expressions ( 20 1 and ( 22 ) in 
p6|. 



9 JO 
irh^^fi^P /A2 

9 W 



dy ■ yjiq 

4q2lnA ) , 



(44) 



as at small y the integral vanishes in both cases. As a result, adding a vortex means adding an amount of 
the Euclidean action 



A^(,).l^^^i!^lnA 



(45) 



8 



Suppose Kq is the fluctuation factor calculated from the path integral in (42) around the g = config- 
uration, and KqKi is that around the q — 1 configuration. If '*"q^'s are the eigenvalues of Ljq, then the 
partition function of the q ~ case is given by |51j 



1 



(46) 



For (7=1, the translation invariance of the vortices gives rise to the existence of the zero modes [52] ■ We 
know that we can generate solutions with w = by simply moving the vortex around. Since the vortex is 
already rotation invariant, it is enough to consider a vortex centered at a; = i?. The displaced vortex solution 
is given by n = ti^ (y'), C = qip' , where 



y 



+ i?2 - 2yR 



For small R we have 



y 



cos if, y snap =y8iiiiip, 



R 



y — R cos Lp,ip' ^ ip -\ sin ip 



(47) 



(48) 



and the deviation from the centered vortex is 

T qR . ^ 

C,i = — sm 1^,0 — —R — ^ cos Lp 

y nq dy 

Then the zero-mode action Sq is given by 

If Aq is the operator for q = and Ai for q = 1, then the fluctuation factor is given by 



(49) 



(50) 



(51) 



where the second equality is due to the existence of a zero mode because of the translational invariance of 
the vortices [HI |32] ) and det' is the determinant excluding the zero mode. The ratio of the determinants is 
given by the Gelfand-Yaglom theorem ^531 [S3] . 

Now consider the situation where more than one vortex is formed. In the dilute gas approximation the 
vortices are assumed to be far apart so adding n vortices of vorticity q — 1 increases the Euclidean action 
by nAT{l) [301 HI]- Form a statistical ensemble of different numbers of vortices n, the partition function is 
given by 



det Ao 


' ISo 


det Aq 


det Ai 


V 2Trh 


det' Ai 



oo „ „ 



y2 . 



d y„ 



Ko exp - 



^o(O) 



A^Kie- 



(52) 



where the integrations over the space have an upper cutoff A^, and the factor ^ is due to the indistin- 
guishability of the vortices. The decay probability per unit time of the configuration from g = to 1 is 

m 



T = -hmT = — Im(ifi)e-^+2i"A. 

h hp 



(53) 



From the expression of the decay probability, the BKT transition temperature can be read off from the 

-jf +2inA ^ ^ given by 



exponential factor since the formation occurs at a reasonable rate as e 



BKT 



2mkf 



(54) 
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where the definition of healing length ^ in Q is used and po = l^/g ,18 is the number density of the lowest 
macroscopically occupied state of the homogeneous configuration — 1. This agrees with the known results 
in the original BKT theory [5] . [56j The correction due to non-homogeneneous configuration in the transition 
temperature is given in Ref . [371 ES] . 

Because Ki is given as the square root of the ratio of the determinants of two differential operators, we 
expect it is of order 1. Then by dimensional analysis, T ^ 7^ ^ 13.1 ms~^. j57| The average time of vortex 
formation is then of the order of 0.08 ms. The numerical estimation of the vortex nucleation rate around the 
transition temperature with the estimated numerical parameters listed in Sec. |ll]is plotted as shown in Fig. 
[2j The vortex formation is very slow below the transition temperature but it increases drastically when the 
temperature increases past the critical point. 



49 - 



30 - 



10 - 



10 - 




FIG. 2: The rate of vortex nucleation around the BKT transition temperature. 



VI. COMPUTING THE IMAGINARY PART OF THE FREE ENERGY 



The expression for the decay rate in ( |53[ ) shows that under the dilute gas approximation the stability of the 
canonical equilibrium hinges on whether the path integral over fluctuations around a one- vortex configuration 
is complex. 



(35). We perform the Gaussian path integration over to obtain 

M/3 



Recall that the action for a linearized fluctuation is given by ( 39 1 , where the Ljq operators are defined in 

in 



S{q) ~ M<l) + 



rfJ-P n 

Jo ''J'yyT. 



j = -oo 



LA 



1 



^3 + g^^jq D'>,-3q^-3)iDs,jq5j) 



(55) 



If the Ljq operators are positive definite, it is clear that the path of steepest descent away from the stationary 
point corresponds to real 5j^ and the path integral is real. 

This is indeed so when we are considering fluctuations around a homogeneous configuration, namely 
q ~ Q, Uq = \. In this case, the eigenvectors of Ljo ^-re Bessel functions of order j. The requirement that the 
Euclidean action must be finite means that we only need to consider eigenfunctions which do not diverge at 
infinity and are regular at the origin. The only Bessel functions satisfying these conditions are of the form 
Jj [ny] corresponding to a positive eigenvalue k^. Thus we conclude that the no- vortex state is stable at 
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zero temperature, when the no-vortex configuration is the only finite action extremal point in the partition 
function. 

Let us see if this argument carries over for nonzero q. For simplicity, we set q = \ (however, we shall leave 
q explicit). We seek finite action solutions to the equation 



h,F^, (y) = -^F^^ (y) , 
with real n. The further change of variables 



(56) 



(57) 



reduces the left hand side to a Schrodinger operator 



where 



V,, {y) 



/« (y) = -n'f^^ (y) , 



(58) 



V,, (y) 



y^ 



a-n, (y)) 



(59) 



Therefore the question of whether the Euclidean action for linearized fluctuations around an isolated vortex 



is positive definite becomes whether a one-dimensional particle of mass 1 /2 in the potential ( 59 1 admits a 



negative energy state. Now, the potential happens to be everywhere positive for all j > 0, so we may discard 
this possibility outright unless j — (see Fig. |3]) 



V[y] 




FIG. 3: The effective potential Eq. ([59} for j 
in Fig [T] to compute Uq. We see that for j 



1 (upper curve) and j = (lower curve). We have used the profile 
1 there can be no negative energy state, even more so for larger values 



of j. For j = 0, on the other hand, the potential is negative for large enough values of y. The existence of a negative 
energy states is shown by Bohr-Sommerfeld condition as in (601. 



In the J = case there is a well defined potential well, and we must investigate whether it is deep enough 
to support a bound state. One possibility is to check the Bohr-Sommerfeld condition, namely, whether there 
is a value of k such that 
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where y± are the classical turning points, namely the roots of k? + Vqi (y) ~ 0. The answer turns out 
to be yes though just barely. Under the approximation given in Fig. [l] for the density profile, the Bohr- 
Sommerfeld condition is satisfied for k = 0.024. The turning points are located at ?/_ = 1.32 and y+ — 21.2. 
Bohr-Sommerfeld quantization would also predict bound excited states; however, these states fall beneath 
the accuracy of our approximations, and they may be considered artifacts. For example, according to Bohr- 
Sommerfeld quantization the first excited state appears at k = 4 x 10^^, with the outer turning point at 
y = 12, 500. This is beyond the intended size of the original homogeneous patch, because from the numerical 

estimation in Sec. |ll| the size of the patch in the dimensionless unit is ^J^^ ^ 40, which is far less than 

12,500. p] 

Observe that not only have we shown that the Euclidean action for axially symmetric perturbations of 
the isolated vortex is not positive definite, but we have also characterized the eigenvector corresponding to 
the direction in configuration space where it becomes negative. Since Uq does not commute with Lqij this 
eigenvector does not correspond to an actual solution of the linearized fluctuations. However, its existence 
is enough to show that the free energy acquires an imaginary part. 



VII. SUMMARY AND DISCUSSIONS 



In this work we have calculated the rate of decay of an effectively homogeneous 2D Bose gas (described 
by T^(x) =0), in the form of Ae~T ^ which complies with the well-known Arrhenius law. The prefactor A is 
proportional to the imaginary part of the fluctuation factor of the free energy of a one-vortex conflguration 
in the path integral. It is known that this imaginary part is due to the negative eigenvalue of the fluctuation 
operator belonging to the eigenvector that defines the direction the fluctuation spontaneously grows along 
(in real time). The qualitative features are like those in the decay of a metastable state due to classical 
fluctuations |38j and barrier penetration due to quantum fluctuations around the instanton solution of the 
Euclidean action [40' . We flnd that the imaginary part comes from the axially-symmetric modes for nonzero 
vorticity conflgurations. As a result, we conclude that while at T = 0, the gas without any vortex is stable, 
the canonical ensemble of different numbers of vortices of the gas is unstable at any flnite temperature. 

Using the fact that at the BKT transition Ae~T ~ 1 we derived the BKT transition temperature Tbkt 
in terms of the number density of the homogeneous phase given in ( 54 1 . This expression derived via thermal 
field theory provides a more quantitative alternative to that originally derived from thermodynamics con- 
siderations of the competition between the energy and the entropy of a vortex [5] . It is known that isolated 
free vortices are formed in the normal phase above the BKT temperature. Hence the decay rate calculated 
here is also the rate of vortex formation. Our calculations show how it increases with temperature. The 
probability of the creation of vortex pairs in a trapped gas increases with temperature as well, as indicated 
by simulation studies with the PGPE [24]." 

In the experiments, measurements on the gas are made some time after the confining potential is abruptly 
turned off. In Dalibard group's experiment there are more isolated free vortices (measured by the dislocation 
of interference pattern of two planes of gas) at higher temperature after 20 ms time of flight (TOF) j32) . 
Our results are consistent with this finding in that the rate of isolated vortex formation increases with 
temperature. In Phillips group's experiment |34| . they observed different characteristics on the density 
profile below and above the BKT temperature after 10 ms TOF, which is at a rate slower than the rate 
of the formation of isolated vortices. However, since we have assumed a homogeneous, time-independent 
configuration as starting point, this should factor in the comparison of our results with experiments. Further 
studies to bridge these gaps are desirable. 
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